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1. Introduction 

The Vlasov equation, often called CoUisionless Boltzmann equation, describes, over 
a certain time frame, the large scale dynamics of Hamiltonian systems of interacting 
particles, in the limit where each particle feels the effect of many others. It is thus 
found in many fields, among which plasma physics of course, and astrophysics, where it 
describes self-gravitating systems. 

The Vlasov equation shows a notoriously rich dynamics. First, it possesses a 
continuous infinity of stationary states, whose study may in itself be a complicated 
problem. The next problem, the study of the linearized dynamics close to a stationary 
state, has a long story. In 1946, Landau [1] formally showed that close to a stable 
homogeneous stationary state, a density perturbation may decay exponentially, which 
was a very surprising result for a Hamiltonian system. This was the starting point of an 
extremely abundant physical literature. On the mathematical side, it has been proved 
that Landau analysis is correct [21 [3] . It is also known that the exponential damping fails 
when the reference state or the perturbation is not analytic |4], or the spatial domain 
unbounded [5]. 

The non-linear case, that is the study of a perturbation close to a homogeneous 
stationary state, under the full Vlasov dynamics, is much more complicated. The subject 
has witnessed spectacular progresses recently PITIIH]- 

Most studies on Landau damping deal with homogeneous stationary states. By 
comparison, the literature on perturbations of inhomogeneous states is less developed, 
although these states are very important. Despite the technical difficulties involved, 
a lot of work has been done in the astrophysical literature. Kalnajs [9], and then 
Polyachenko and Schuckman [10] have developed a powerful formalism, sometimes 
called the "matrix method", to solve the linearized Vlasov equation in a inhomogeneous 
context. Since then, it has been used to study the instability of many models of stellar 
systems, and compute the growth rates (see [11] for a textbook account). Some other 
methods to compute instability rates have been introduced and applied to toy models 
recently [T2l [T3l [T4] . Purely oscillating modes in ID have been investigated in [15]. Also, 
when the question is to prove stability and not to compute a decay rate, there exist 
powerful variational methods [H]. However, computing decay rates for non-oscillating 
stable situations is more tricky than computing instability rates, since it involves an 
additional analytic continuation, as in Landau's original analysis. In the astrophysical 
literature, we are aware of only two papers performing this continuation numerically, 
and thus explicitly computing the analog of Landau damping rates [16], [17]. In [TH] . 
the continuation is performed analytically, but on a one dimensional toy model. Thus, 
although Landau damping is considered an important player in the dynamics of stellar 
systems (see for instance [H]), there are very few actual computations of damping rates 
close to inhomogeneous stationary states. 

In addition, these studies done in the astrophysical context do not mention a 
fundamental difference between homogeneous and inhomogeneous stationary states: in 
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the inhomogeneous case, no matter what is the regularity of the stationary state and 
the perturbation, the asymptotic hnear decay is never exponential. This has been seen 
in the Kuramoto model, which shares some properties with the Vlasov equation [T9] . 
and on a type of Vlasov equation in [20]. The picture is as follows: the dynamics may 
show a transient exponential decay governed by a Landau pole, but the asymptotic 
decay is always algebraic. This phenomenology is also well known for the 2D Euler 
equation [21], which is in many respects similar to the Vlasov equation. We note that 
such algebraic decays may also arise close to a homogeneous state, when the perturbation 
or the reference state is not regular enough [1]; we stress that the situation is different 
close to inhomogeneous states, because the algebraic decay occurs for all reference states 
and perturbations. [22] contains a detailed analysis of this phenomenon in ID, including 
derivation of the decay exponents and comparison with direct numerical simulations. 
In this paper, we extend the analysis of [22] to multidimensional systems, putting our 
emphasis on 2D systems and 3D systems with spherical symmetry. These includes some 
important models of stellar systems. 

More specifically, we first use the standard matrix method to formally solve the 
linear dynamics in a Laplace transformed space (section |2]). From this starting point: 

i) we identify and classify the singularities appearing in the Laplace transform of the 
perturbation. It turns out that the zoology of singularities in two dimensions is much 
richer than in one dimension, see section |3l 

ii) we exhibit the asymptotic decay for each type of singularities in section |H 

iii) we introduce advection equations associated with the linearized Vlasov equations 
in section [51 These advection equations have a solution in integral form; hence the 
numerical task to obtain the temporal evolution of the system is reduced from solving 
the (linearized) Vlasov equation in 2(i-dimensional phase space to performing integrals 
in a dimensional space, where d is the spatial dimension. 

iv) Using the above theory, we analyze the asymptotic decay of a perturbation in a 
toy model and a spherically symmetric stellar model, and check the theory against the 
numerically computed exact solutions to the simple advection models, see sections 16.41 171 

Our analysis of the asymptotic decay of a perturbation is formal, and moreover 
relies on the linearized Vlasov equation. There is no guarantee that this is relevant to 
understand the asymptotic decay of the non-linear Vlasov equation. In principle, our 
results should then be supplemented by direct numerical simulations of the full Vlasov 
equation. However, such simulations are difficult in more than one spatial dimension, 
since we are aiming at an asymptotic in time regime, while keeping a good spatial 
precision. We left this study for future work. This is why we analyze instead with the 
present theory easily solvable linear advection equations, for illustrative purpose. 

2. Solution to the linearized Vlasov equation: the matrix method 

Throughout the paper, we will use bold letters (q, p, 6, J, . . .) to denote vectors. 
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2.1. System 

We start with the Vlasov equation in d spatial dimensions for the one-particle 
distribution function f{q,p,t), 

dj + VpH ■ V J - V,H ■ Vpf = 0, (1) 

where qr G X C M'^ is the position variable, p G M"' the conjugate momentum variable, 
and Vq, Vp denote respectively the gradients with respect to q and p. The domain X 
is in a stellar system, and [0, 2'irY if the system has periodic boundary conditions. 
H is the one-particle Hamiltonian defined by 

H[f] (q, P,t) = ^ + t) + ^^q). (2) 

The potential is defined from the two-body interaction potential v and the 

distribution / as 

$[/](q,t)= [ dp [ dqv{q-q')fiq',p,t). (3) 
The external potential $ext(9) creates a force F(,y±[q): 

FoAq) = -Vq$ext(g). (4) 

Let /o be a stationary solution to the Vlasov equation ([1]). For the stationary 
solution, the potential $[/o] and the one-particle Hamiltonian H[fo] are independent of 
time. In this paper, we limit ourselves to situations where the one-particle Hamiltonian 
H[fo] is integrable. Thus, we may introduce actions J = (Ji,..., J^) and conjugate 
angles 6 = {6i, . . . ,6^), and the one-particle Hamiltonian is a function of the actions 
only, i.e. iJ[/o](J). A practically important example is given by spherically symmetric 
stellar systems, see section 13.11 A stationary solution can be constructed by taking /o 
as a function of actions, satisfying self-consistent conditions, since the actions depend 
on /o through the potential $[/o]. 

We now consider a small perturbation to a stationary solution fo{J): 

f{e,j,t) = fo{J) + fi{e,j,t) (5) 

and write the linearized Vlasov equation for /i: 

dJ, + n{J) ■ Veh - Vjfo ■ Ve$i = , (6) 
where $1 = is the perturbed potential 

^iiq,t)= [ [ viq-q')Mq\p,t)dq'dp, (7) 

and fl{J) = Vji/[/o](J) is the vector of the frequencies in the unperturbed potential. 
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2.2. Fourier- Laplace transform 

To analyze the linearized Vlasov equation ([6]), we introduce the Fourier-Laplace 
transform ^(m, J, uo) of a function J, t) as 

/TT /' + 00 

rf^ e^^"*-^ / dt e''^'u{e,J,t) (8) 
TT Jo 

where m = (mi, . . . , m^) is a d-uplet of integers, m ■ 6 is the Euclidian inner product 
and Im(a;) is large enough to ensure convergence of the integral with respect to t. The 
inverse transform is then 

n{e, J, t) = — ^ Yl I ^(^' (9) 

where F is a Bromwich contour running from — oo + ia to +oo + io", and the real value 
(T is larger than the imaginary part of any singularity of u{m, J,u) in the complex 
cu-plane. 

Simple and standard algebraic manipulations on ([6]) yield the following expression 
for the Fourier-Laplace transform of the perturbation: 

/i(m, J, uj) = A{m, J, co')<l>i(m, J, u) + B{m, J, cj), (10) 

where 

-^^"^"''"^ = m.fl(J)-c. ' 

S(m,J,a;) = ^^^^^, (12) 

and ig{m,J) is the Fourier transform of the initial perturbation fi{0,J,t = 0) with 
respect to 6. Without loss of generality, we assume 

(7(0, J) = 0. (13) 

Notice that we can always include a non-zero g{0, J) in the stationary state /q. 

2. 3. Biorthogonal functions 

We should now solve the two coupled equations ([7]) and ffTOj) to find /i and $i. The usual 
strategy at this point is to introduce two families of biorthogonal functions O [101 [23] ■ 
The computations below are standard (for a textbook reference, see for instance |llj): 
we reproduce them here to ensure that this article is self-contained. 
We expand any density function p{q) on the basis {di{q)}i<^j: 

p{q) = Yaidi{q) (14) 

i&I 

and any potential function $(g) on the basis {uk{q)}k&K'- 

Hq) = J2bkUkiq). (15) 

k&K 

The two index sets / and K satisfy K G I G 1j, and the two basis {c?j}jg/ and {uk}keK 
are chosen such that 
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(i) the two families are orthogonal to each other: 

{di.Uk) = / di{q)uk{q)dq = \Ak, e /, k ^ K) (16) 

with Afc 7^ 0, and where 5ik is the Kronecker 5. 

(ii) Mfc is the potential created by the density distribution di. 

Assuming we have at hand two such families of biorthogonal functions, we expand 
the perturbation density 



pM.t) = / h{q,p,t)dp (18) 
and potential as 

pi{q,t) = ^a,(t)c?,(g) (19) 

$i(g,t) = ^afc(t)Mfc(q). (20) 

k&K 

Properties f|T6l) and f|T71) ensure that the coefficients are identical for pi and $i. 
Using (l20l) . we first compute $i, the Fourier-Laplace transform of $i: 

<li(m, J,a;) = ^afc(a;)cfc(m, J) (21) 

where dk is the Laplace transform of a^, and Ck{m, J) is the Fourier transform of m^: 

Cfc(m,J)= / Uk{q)e-'^-''de. (22) 



Substituting fl^T]) into fITU]) . and performing the inverse Fourier transform, we obtain 
an expression for fi{6, J^u), the Laplace transform of the perturbation: 

Now, we multiply ( l23l) by ui{q) (/ G K) and integrate over dOdJ. We compute the left 
hand side using the change of variables {6, J) — > (q,p) which has Jacobian 1, and the 
property (flGl) : the result is \idi{uj). In the right hand side, performing the integration 
over 6 introduces the functions ci{m, J). The result is 

Xidiiio) = J2 Fik{io)dk{uj) + Gi{uj) (24) 



k&K 



with 



and 
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The contributions from m = clearly vanish for F^^, and also for G^, thanks to 
assumption f|T3|) . Defining the {^K) x {^K) matrices A and F{uj) = {Fii;{uj)), where 
A is diagonal with elements {\i}i^k, and the (tl/'^)-dimensional vectors G{uj) = {Gi{u)) 
and d{uj) = {ak{u)), the above equations may be rewritten in compact form: 



The equation det(A — F{u))) = is sometimes called the dispersion relation. 
2.4- Analytic continuation 

Assuming that the decay for large J is fast enough in the integrands, and that the 
Ck{m, J) are regular, expressions (l25l) and (l26l) show that functions Fik{uj) and Gkioj) are 
analytic in the upper half plane Im(a;) > 0, since the corresponding integrals over J have 
no singularity. The integrands are singular however for any real uj at which m-Q,{J)—uj 
vanishes. Thus, expressions ( l25l) and ( l26l) should be analytically continued to define the 
F's, G's and a's in the lower half plane Im(u;) < 0. This analytical continuation is a 
generalization of the usual Landau prescription in ID. It is sketched in some details in 
the 2D case in section 13.21 

The last step to compute the evolution of the perturbation is to perform an inverse 
Laplace transform on the functions a{u:). The large time behaviour of ak{t) is determined 
by the singularities of akioj): our goal now is then to study and classify the singularities 
of a{uj). 

2.5. Singularities of a{uj) and roots of the dispersion relation 

The singularities of a may come (i) from the roots of the dispersion relation det(A — 
F{(jj)) = 0, and (ii) from singularities of the functions F's and G's themselves. We 
discuss the former singularities in this section: we will see that these singularities do 
not dominate in the asymptotic regime. The dominating latter singularities are classified 
in the next section in the two-dimensional setting. 

Roots of the dispersion relation det(A — F{uj)) = yield poles for the functions 
dj. Such poles in the upper half plane correspond to eigenvalues of the linearized 
Vlasov operator, with exponentially growing eigenmodes. Since we are interested in the 
relaxation of a perturbation close to a stable stationary state of the Vlasov equation, we 
assume that there are no such eigenmodes. Poles on the real axis correspond to purely 
oscillating eigenmodes. In order to study the decay of perturbations, we also assume 
that there are no such modes. There may be roots of the dispersion relation in the lower 
half plane Im(ci;) < (in this case, they are rather roots of the analytic continuation of 
the dispersion relation). These are the usual "Landau poles", giving rise to exponential 
damping. This damping may be an important feature of the dynamics at intermediate 




(27) 



The formal solution for a is 




(28) 
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time scales, especially if the pole is close to the real axis [121 HZJ dS], but the asymptotic 
regime is always dominated by the singularities of the functions F's and G's, as will 
become clear in the following sections. 

3. Singularities of F and G in two-dimensional setting 

We have seen that the functions F's and G"s have no singularities for lm{uj) > 0. The 
integrands of F's and G"s are singular for J s such that m ■ n{J) — u vanishes, and 
these singularities may yield branch points on the real u axis. Thus, to investigate 
the asymptotic behaviour of a perturbation due to these branch points, we turn now 
to the singularities on the real axis. From now on and for simplicity, we restrict 
ourselves to two-dimensional integrations over J = ( Ji, J2) in the formulas fl25|l and fl26l) 
entering into the definitions of F and G. We have in mind perturbations in spherically 
symmetric self-gravitating systems, which fit into this two-dimensional setting, as shown 
in section 13. 1[ There would be no major obstacles to an analysis in higher dimension, 
except the growing number and complexity of the possible types of singularities. 

3.1. Spherically symmetric system 

Although Vlasov-type equations appear in various settings, the main situation we have 
in mind is the dynamics of a perturbation of a spherical self-gravitating system. The 
dynamics in the potential created by the stationary distribution is then integrable. It 
is customary [H] to use as actions Jr ("radial action"), L (modulus of the angular 
momentum) and (projection on the z axis of the angular momentum). The radial 
action Jr, energy of a particle and angular momentum L are related by 



where $(r) is the potential created by the stationary state and r^m and r^ax are the 
minimum and maximum distance from the origin reached by the particle. In this 
coordinate system, the one-particle Hamiltonian depends only on the actions Jr and 
L. 

Let us assume that /o and g depend on the actions only through the one-particle 
Hamiltonian. Then, fo{J), g{J) and Ck{m,J) depend only on actions Jr and L, and 
hence the potentially singular integrands in fl25|) and fl26l) also depend only on Jr and 
L. The three-dimensional integrals are thus reduced to two-dimensional integrals. The 
case of spherically symmetric systems is then included in the abstract setting introduced 
in section 13.21 




(29) 
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3.2. An abstract problem 

The abstract problem is to study the singularities of the analytic continuation of <^{z) 
defined for Im(2;) > as integrals over a domain D dM?: 



with fi and ly real functions. Functions F and G, 0251) and fl26l) respectively, fit in this 
framework, with fi{J) = m ■ fl{J), and by defining properly We assume m 7^ 

since the contributions from m = in the definitions of F and G vanish. We also 
assume that n and z/ are very regular: although they are naturally defined over the 
domain D C M^, they may be analytically continued over the complex J domain. We 
further assume that their integrability properties are as good as needed. Our goal is to 
study the singularities of the function 

= lim ip{x + iy), with x,?/GM. (31) 

This is the analytic continuation on the real axis of y:>{z), which is a priori defined for 
lm{z) > 0. 

For z real, the denominator in fl30l) . — z, may vanish. We explain 

in Appendix A why for a generic z there exists an analytic continuation of <y9 in a 



neighborhood of z, and hence there is no singularity for 0. From the computation 
in Appendix A special values of xq = corresponding to a singularity for </> are 



easily identified; we classify the special points J* and the associated singularities for (j) 
in the following subsections: 

(i) Vertex singularity: a point J* at which the boundary of D is not regular 
(section 13. 3p 

(ii) Tangent singularity: a point J* at which the level set fi{J) = fJ'{J*) is tangent to 
the boundary of D (section 13.41) . 

(iii) Critical singularity: a critical point J* of the function fi{J) inside D (section l3.5p . 
We will compute singularities of </)(a;) around Xq = ^i{J*) by expanding /i(J) around 

r. 

fi{J) = Xo + lii{Ji - Jl) + li2{J2 - J2) + ^f^niJi - Jl? 

+ \^^22{J2 - j;? + ^^l2{■h - JW2 - ^2) + ■ ■ ■ (32) 

where, for instance, yUj (z = 1,2) represents {dfi/dJi){J*) and is not zero in general. 
The above three types of singularities (i), (ii) and (iii) appear at the point J* where 
respectively (i) neither fii nor /i2 vanish, (ii) /ii or ^2 vanishes, but not both, and (iii) 
both fii and ^2 vanish. The three types of singularities have therefore respectively 0, 1 
and 2 analytic conditions for fi{J), and 2, 1 and geometric conditions; see Table [TJ 
Considering that there are two variables Ji and J2, these singularities appear generically. 
Special treatments are devoted to the case of an infinite domain D (section 13.61) . and 
to a special situation which however does arise generically for spherically symmetric 
self-gravitating systems (section [3?7I) . 
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Singularity 


Conditions for /u(J) 


Geometric conditions at J* 


Vertex 




On dD, non- regular point 


Tangent 


/xi = or /i2 = 


On dD 


Critical 


/ii = and /i2 = 





Table 1. Summary for the vertex, tangent and critical singularities, /i^ represents 
{df-L/ dJi){J*) for i = 1,2, and is not zero in general. The symbol dD denotes the 
boundary of D. The "or" in tangent singularity is exclusive. 



3.3. Vertex singularity 

We consider first a singular point J* on the boundary of the domain Z), such that there 
is a jump in the boundary's slope at J*. We refer to such a point as a "vertex", see 
figure [TJ This type of singularity is always present in spherical self-gravitating systems 
at zero radial action and zero angular momentum. 

To investigate the leading singularity of 0(x) created by the vertex, it is enough to 
keep only the leading order terms in the expansions of and We expand the 

function yu(J) as 

/i(J) ~ Xo + /il(Jl - J*) + /i2(^2 - ^2*) (33) 

where Xq = fi{J*), and iy{J) as 

iyiJ)r^CiJi-J*riJ2-j;r\ (34) 

where ai and 02 are non-negative integers. We assume that the first-order derivatives 
fii and fi2 do not vanish, which should be the generic case. Here and in the following C 
is a constant whose value plays no role and may vary from line to line. 

The singularity of 0(x) is then investigated on the following reduced expression for 

where Uj* is a domain including a vertex at J = J*. Shifting and rescaling Ji and J2, 
we have 

ip{z) = C / -d.hdJ2, (36) 



juo Ji + J2- (z- Xo) 
where Uq is a domain including the shifted vertex at J = 0. Using new variables 
u = Ji + J2 and V = Ji — J2, 



h I 



(p{z) = y Cki / -, -dudv, (37) 

JtT U — [Z — Xn) 

k,l>0,k+l=ai+a2 ^ ' 

where Cki is a constant and [/ is a domain including a vertex at {u,v) = (0,0). We 
assume that the line u = does not coincide with the boundary of U, since this case 
results in the line singularity (see section [3?7|) . As shown in Appendix B, the integral 
over V yields 

^{z) = C / -du, (38) 

Jo u-{z- Xo) 
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Singularity at infinity 

--,,7 



Figure 1. Vertex (top left), tangent (top right), critical (bottom left) singularities, 
and singularity at infinity (bottom right). In all cases, the domain D is shaded, and 
the dashed lines are level sets of the function /-t(J). 




which fits in the framework (lE.2p . Therefore, the singularity of at x = xq is: 



4f^l^{x) = Ci(x - Xo)°^+"^+^ In |x - xol + C2{x - Xof'^''^+^H{x - Xq), (39) 
where H is the Heaviside step function. 

3.4- Tangent singularity 

We consider now a regular point J* on the boundary of the domain D such that the 
level set = = is tangent to the boundary of D at J = J*; see figure [H 

The appearance of a tangent singularity requires that one derivative of /i vanishes, at 
a point located on the boundary of D. Changing variables, we may assume that the 
boundary of D is locally parametrized as J2 = J2', then the tangency condition reads 

dfj. 



= 0, J* G dD (40) 



dJi 

and at leading order we will use 

/i(J) ~ xo + ^/iii(Ji - J*if + /i2( - ^2)- (41) 

Using also the leading order approximation ( !34l) for i>, and introducing these into ( l30l) . 
we obtain after shifting and scaling Ji and J2 the reduced expression for f{z): 

p jai ja2 
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where the new tangent point is J* = (0, 0), and we have chosen a subdomain U close to 
J* . Changing variables again and performing one integration, the function ip is further 
reduced to 

rc l/2+ai/2+a2 

^{z) = C{1 - i-iy+n / 7 r^^, (43) 



see 



Appendix C for details. This form fits in the framework (]E.2P : hence the singularity 



of (h is: 



Ci{x - Xo)^/2+''i/2+''^i^(x - Xo) 

ai. odd. 



3.5. Critical singularity 



We consider now a point J* belonging to the interior of D such that the level set 
/i(J) = fi{J*) = Xq is singular at J*; this requires that both derivatives of /i vanish at 
J*. In addition, we assume that this critical point is generic, thus the Hessian at J* has 
both eigenvalues different from 0. We also discard as non-generic the situation where the 
critical point J* belongs to the boundary of D; see an illustration on figured] The critical 



point may be a local extremum or a saddle point. We study both cases in [Appendix D 
and summarize the results here. Using the leading order approximation f p^ for u, 
changing variables and performing one integration, the function ip{z) is reduced to 

rc (ai+a2)/2 

^{z) = C / -du (45) 

Jq U-[z- Xq) 

for a local extremum, and to 

/c „ (ai+a2)/2 1„ L.l 
, \^ du (46) 
U-[Z- Xq) 

for a saddle point. We note that these expressions are valid for even ai and 02, and 
no singularity appears otherwise. These functions fit into (IE. 21) and (IE. 31) respectively, 
and, for both extremum and saddle points, the singularity of is: 

Ci{x - xo)("i+"2)/2 1^ 1^ _ 
/,Jg(2;) = <j +c'2(x - xo)(''^+'^^)/2i/(x - Xo) ai, as : even (47) 
otherwise. 

3.6. Singularity "at infinity" 

Assume yu(J) tends to when \ J\ tends to infinity; then xo = is a singular point for 
0(x), see figure [H This is a common situation for spherical self-gravitating systems. To 
study this singularity, we assume that the domain D is [0, +00) x [0, +00), and 

^l{J) ~ (Ji + J^)-"" , v{J) ~ (Ji + J2)~'. (48) 

We assume that both a and h are integers with 6 > 2 to ensure convergence of the 
integrals in ( l30l) . 
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Changing variables to u = Ji + J2, v = Ji — J2, and then from u to s = d^/u^, we 
obtain 

pc (b-2)/a-l 

= Cl -ds. (50) 

Jo S-[Z- Xq) 

The integral over s fits in the framework (lE.2p . and the leading singularity of at Xq is: 

' Ci{x - XoY^''^^/''~^\n\x - Xo\ 

+ C2{x-XoY''-^^/''-^H{x-Xo) 6-2 = (mod a) 
Ciix-XoY'-^^/''-^H{x-Xo) ^ ^ 

+ C2{xo - x)(^-2)/»-ii7(xo -x) (mod a). 

We will find this kind of singularity in a model of self-gravitating systems, see section [7l 

3.7. Line singularity 

We consider now a case which seems at first sight very special, but which does occur 
generically in spherical self-gravitating systems. In a spherical potential, the orbit of a 
particle with angular momentum L = is purely radial; for small L, it is very elongated, 
and remains close to be purely radial. When the particle travels from one maximal 
radius to the next along the elongated orbit, this corresponds to one radial period; this 
also corresponds approximately to a half angular period, see figure O Thus, on the 
L = line, the radial frequency is exactly twice the angular frequency fir = 20,0] or 
equivalently m ■ fl = 0, with m = (1, —2). 

To study this singularity, we consider a domain D as in figure E) the axis Ji = is 
a boundary. Suppose that the function yu( J) is constant on the J2 axis as 

/i(0, J2) = const. = xo , VJ2 > 0. (52) 

We then expect a singularity of 0(x) a.t x = xq. Close to the singular line Ji = 0, we 
expand /i(J): 

l2{J)=Xo + W{J2)Jl + ... (53) 

Always having in mind self-gravitating systems, we assume that it is possible to expand 
u close to the J2 axis as: 

z/(J) ~ /i(J2)Jr (54) 
with a non-negative integer Oi. The function (p around z = Xq reads 

= rdJ, r^J2 „,,,y^f r (55) 

Jo Jo W{J2)Ji-{z-xo) 

The integral over u fits in the framework (1E.2P : we conclude that the singularity of at 
Xq is: 

0^°g(a;) = Ci{x - XoT' In \x - Xo\ + C2{x - x^Y^H{x - Xq). (57) 
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Figure 2. Left: schematic picture of a line singularity. Right: orbits of a particle in 
a spherical potential with L = (black line) and L small (blue elongated orbit). 



This computation relies on the fact that W{J2) is well behaved (for instance bounded 
away from and oc). 

4. Asymptotic behaviour of a perturbation 

The singularities in the functions F and G are passed on to the coefficients dk{u}). Each 
singularity corresponds to a decaying term in the inverse Laplace transform. If ak{oo) 
had a finite number of singularities, the asymptotic in time behaviour of ak{t) would be 
a sum of decaying terms, each one corresponding to a singularity of dk{oj) (see Theorem 
19 in [21]). We expect that this remains true in our situation, where dk{uj) actually has 
an infinite number of singularities. 

According to the analysis of section [3] the singular part of dk{u}) (or rather of its 
analytic continuation on the real axis) close to a singular value Uq is, for u > Uq: 

{u — ojq)^ In |a; — Wol for n G Z+, 
or (w - ujQfH{uj - Wo) for a G M \ Z. 

Here Z_|_ is the set of non- negative integers. The behaviour is of course similar for u < uq. 
If a function (p presents any of the above singularities, we associate to it an asymptotic 
behaviour of its inverse Laplace transform 0, according to the following rules [2^ : 

CfH = ^(^-^orin|a;-c^o| (r^GZ+) ^ 0-(t) = , (58) 

Cr(^) = Cicu - UorH{u -Uo){aeR\Z)^ 0-(t) = C'^ . (59) 

Singularities and associated asymptotic behaviour are summarized in Table [21 Thus, to 
extract the asymptotic behaviour of afe(t), one needs to 

(i) Enumerate the singularities appearing in dk{uj) 

(ii) Keep the strongest singularity, corresponding to the slowest decay 
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Check if the symmetries of the problem induce a special cancellation. 

Unless a cancellation occurs, this strongest singularity yields the decay exponent and 
asymptotic frequency of akif). We illustrate this strategy on several examples in the 
following sections, including the case of a cancellation; a careful reexamination of the 
singularity is then needed. 



Type 


Singularity 


Exponent a 


Damping 


Sign 


Vertex 


flE.4p 


1 + oi + a2 


g— i2;oi^-(2+ai+a2) 




Tangent 


(ai : even) 


(|E.5P 


l/2 + ai/2 + a2 


g— i2;o<^-(3/2+ai/2+a2) 






(ai : odd) 


N/A 


N/A 


N/A 




Critical 












extremum 


(ai, 02 : even) 


(|E.4) 


(ai + a2)/2 


g-ixoi^-(l+(ai+a2)/2) 




saddle 


(oi, a2 : even) 


(]E.6|) 


(ai + a2)/2 


g-ia:oi^-(l+(ai+a2)/2) 


(-1)" 




(otherwise) 


N/A 


N/A 


N/A 




Infinity 


((6-2)/a e Z) 


(IE.4P 


(6- 2)/a- 1 


^-ixotj.-(b-2)/a 






((6-2)/a^Z) 


(1K5) 


(6-2)/a- 1 


g-ixoi^-(b-2)/a 




Line 


flE.4|) 


ai 


g-ixot^-(l+ai) 





Table 2. Table of singularities of the function (j){x) defined by pO|) and ((3T|) . with 
their associated exponent a (see section[3land |Appcndix E ), and asymptotic dampings. 
The singularity is produced at J* = {Ji, J2), and xq = /.((J*, Jj). The leading order 
of the numerator is v{J) — {Ji — Ji)°'^{J2 — J2)"'^ for vertex, tangent and critical 
singularities, I'iJ) — ( Ji + ^2)^^ for a singularity at infinity, and i^(J) ~ /i(J2) J"^ for 
a line singularity. We assume /i(J) — ( Ji + ./2) for a singularity at infinity, ai and 
02 are non- negative integers, a is an integer and b is an integer satisfying b > 2. N/A 
means that no singularity appears in 4>{x). This table shows the leading singularity 
for each type. As shown in sections |6] and [71 a special cancellation for the leading term 
may happen between modes m and m. The column Sign represents the relative sign 
between the singular parts of modes m and— m for xq = 0. No cancellation is expected 
for Xq irrespective of the relative sign. A bar means no simple relation between 
the two modes. Note that the relative sign may also depend on the numerator i^(«7) of 
([30]) if it depends on m, but this dependence is ignored in this table. 



A remark is in order: we study here the asymptotic behaviour of the ak{t), which 
are the coefficients of the density and potential perturbation in the expansions (|T9l) 
and (I20I) . One might expect that the slowest decay for these coefficients yields the 
asymptotic decay of the density or potential perturbation at a given point in space. 
This is not necessarily true, as we have no information on the rate at which reaches 
its asymptotic regime, and thus have no control on the series (fT9|) and (l20l) . 



5. Preparation for numerical tests 



5.1. On numerical tests 

The theory developed in this paper should describe the asymptotic behaviour of a 
perturbation evolving according to the linearized Vlasov equation. Two questions arise: 
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i) Since the above computations are formal, one may wonder if the theory is correct for 
the evolution of a linearized Vlasov equation. 

ii) Even if the theory is correct for the linearized Vlasov equation, it is not clear that it 
describes the asymptotic behaviour of the fully non-linear Vlasov equation. 

To answer these questions, direct numerical simulations of the Vlasov equation are 
required. Such a check has been done in a one-dimensional setting using iV-body 
simulations [22]; however, the memory, accuracy and time frame required to perform 
similar simulations in two or more spatial dimensions make the task challenging. Instead, 
we will illustrate the various aspects of the theory on much simpler linear advection 
equations in section [6], and apply it to make definite predictions on a spherically 
symmetric self-gravitating system in section [3 

5.2. Advection equation and singularity analysis 

To illustrate the theory, we consider an advection equation 

dtf + n{J) ■ Vef = 0. (60) 
The exact solution is 

f{e,j,t) = f{e-n{j)t,j,t = o). (6i) 

We consider the temporal evolution of the expected value of A{6) defined by 

{A){t) = J AiO)f{e,J,t)dOdJ, (62) 

where / is governed by the advection equation fl60|) and we assumed that A is a function 
of only. From the Fourier-Laplace transform of the advection equation ( l60|) . the 
Laplace transform of (A) (t) is given by 

^ E / me'-'^O I J^(f^\j J, (63) 

where ig{m, J) is the Fourier transform of f{0, J,t = 0). The function 

^(a;;m)=/ j^[r:^'f dJ (64) 
J m ■ il[J) — 00 

fits in the framework fl5Dl) . thus the asymptotic damping of the expected value {A) (t) 
is predicted by the theory presented in sections [3] and HI 

The theoretically predicted asymptotic damping can be checked without computing 
the temporal evolution of / numerically, since we have the exact solution fl6T]) to the 
advection equation, which gives 

{A) it) = J2 [ A{e)e''^-^dO I ig{m,J)e^'"'-"^-^^'dJ. (65) 

m 

Integrals over can be performed analytically, and hence our numerical task is to 
perform integrals over J. Consequently, we can reduce the 2 d- dimensional problem fl^Ul) 
to (i-dimensional integrals over J. For 2-dimensional cases, numerical integrations are 
possible with good accuracy. 
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6. Numerical tests 

We set the initial distribution as 

f{e,J,t = 0) = jfV^(Ji-J*)"n^2-^2)"'e-('^i+''')/2(l+cos^i cos^2), (66) 

where the unity is to keep / positive, but is not important since the mode (0, 0) is 
temporally invariant. We therefore focus on the four modes m = (1, 1), (1, —1), (—1, 1) 
and (—1,-1). We assume that the system is defined on the domain (Ji, J2) ^ [0, +00)^. 
The factor (Ji — J*)"^(J2 — ^2)"^ corresponds to the leading estimation of i^(J) (IMl) . 
and the factor will be used to select the singularities we want to illustrate. The 

function g is the same for the four modes and, remembering ig is the Fourier transform 
of the initial perturbation, 

^g{m,J) = ^J^^J^-(J,-J*r(J,-j;y-e-^'?^-''^^/\ m = (±1,±1).(67) 
We introduce four observables: 

Ai = cos(^i+^2), A2 = sin(^i+^2), ^3 = cos(^i-^2), ^4 = sin(^i-^2), (68) 

to pick up singularities appearing in the modes (1, 1) and (—1, —1) by Ai and A2, and 
in the modes (1, —1) and (—1, 1) by A^ and ^4. 

We now report detailed theoretical predictions and numerical illustrations for 
vertex, tangent and critical singularities in sections 16. ![ 16.21 and 16.31 respectively. 
These include cases with cancellations between two modes. The line singularity is 
discussed in section 16.41 on a toy model, which also has the vertex and tangent 
singularities. We discuss the dominant singularity for different modes, again taking into 
account cancellations, and check the predictions against numerical computations. The 
singularity at infinity will be demonstrated in section [7] on a spherical self-gravitating 
model. 

In this section, all the numerical integrations for the exact expected values are 
performed by introducing a cutoff at Ji = J2 = 10 and slicing the interval [0, 10] in 2^'^ 
bins. 

6.1. Vertex singularity 
We set the frequencies as 

= Ji, n2{J) = J2, (69) 

which yields a vertex singularity at the origin, with (^1,^2) ~ (0,0); the associated 
frequency is uq = m ■ fl{J*) = 0. We compute {Ai) (t) for ai = 0, 1 and 2 with fixed 
02 = 0. There are no singularities of other types, thus we may use /ii = /12 = 0. From 
Table O we predict dampings and cancellations as shown in Table [3l Cancellations 
occur for Ai and A3 with an odd ai, since the cosine observables pick up the sum of the 
modes m and —m; for instance: 

= cos(^i + ^2) =^ (ir)(c^) oc<^(a;; (1,1)) + (^(o;; (-1,-1)), (70) 
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(01,02) 


(0,0) 


(1,0) 


(2,0) 


(^1) (t) 




c 




w it) 


c 


r3 


C 


{A3) it) 


r2 


c 




(A,) it) 


c 


t-' 


c 



Table 3. Dampings and cancellations for the vertex singularity. "C" stands for 
" cancellation" . 




Figure 3. Dampings for the vertex singularity. Observables are (a) Ai and (b) A2. 
The parameter 02 is fixed as 0. Observables ^3 and A4 give qualitatively similar panels 
to Ai and A2 respectively. 



and (fiioj; (—1,-1)) = i—l)°''^^"'^(piLj; (1,1)) holds for the vertex singularity. Similarly, 
cancellations occur for the sine observables A2 and A^ with an even ai since, for instance, 

= sin(0i + ^2) =^ (X)(w) (X(/.(a;;(l,l))-(/.(a;;(-l,-l)). (71) 

We note that the cancellations occur for the leading singularities; higher-order 
singularities may survive in general. However, in this case, higher-order singularities 
coming from the expansion of exp(— (Jf + J|)/2) do not change the parity of Oi + 02, 
since expanded terms consist of Jf and J|. Thus, all orders are cancelled when the 
leading order is cancelled. 

The above prediction is confirmed by direct numerical integration of f l65|) as shown 
in figure El We remark that dampings without oscillation reflect the zero frequency 
Uq = 0. We have also checked that no cancellation occurs when we replace Qi by 
fliiJ) = Ji + 1, since = 1 7^ (figure not shown). 

6.2. Tangent singularity 
We set the frequencies as 

^]l(J) = (Jl-l)^ n2iJ) = J2. (72) 

There are a vertex singularity at ( Ji, J2) = (0, 0) and a tangent singularity at ( Ji, J2) = 
(1,0). We are interested in the tangent singularity, and hence we set (^1,^2) — (1,0) 
which gives the frequency ojq = 0. We choose ai and 02 from the set {0,1,2}, and 
use hi = 2 and /i2 = to hide the vertex singularity: the slowest decay coming from 
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the vertex singularity is then t~^. The predicted dampings and cancellations are shown 
in Table m An explanation is needed for the case (01,02) = (1,0): the leading order 
does not give any singularity, since ai is odd; however, expanding ig{m, J) around 
( Jj*, J2) = (1, 0), the second leading order gives ai = 2; thus the damping must be 
as in the case (01,02) = (2,0). We remark that no cancellation occurs for the tangent 
singularities. These predictions have been checked numerically, though no figures are 
reported. 



(01,02) 


(0,0) 


(0,1) 


(0,2) 


(1,0) 


(2,0) 


(A) (t) 


^-1.5 


^-2.5 




r2-5 


^-2.5 



Table 4. Dampings and cancellations for the tangent singularity. 



6.3. Critical singularity 
We set the frequencies as 

fil(J) = (Jl-l)^ n2iJ) = iJ2-iy. (73) 

Singularities vertex singularity at (Ji, J2) = (0,0), tangent singularities at 

{Ji, J2) = (1,0) and (0,1), and a critical singularity at {Ji, J2) = (1,1)- We are 
now interested in the critical singularity, and hence we set (^i , J|) = (1,1) giving 
the frequency = 0. We choose Oi from the set {0, 1, 2}, and use /ii = /12 = 2 to hide 
the vertex and tangent singularities as done in section 16.21 

Cancellations between the modes m and —m are as follows. We first remark that 
no leading singularity appears for (01,02) = (1,0) since oi is odd. In this case, as 
discussed in the case of the tangent singularity, the second leading singularity gives 
the same singularity as (01,02) = (2,0), thus we concentrate on (01,02) = (0,0) and 
(2,0). The second remark is that the critical singularity at (Ji, J2) = (1,1) behaves 
as an extremum singularity for Ai and A2, and as a saddle singularity for and A4. 
Referring to Table |2] and the above second remark, the leading singularities in Ai and 
are cancelled for (01,02) = (0,0), and in A2 and ^3 for (01,02) = (2,0). In the 
latter case, the next leading singularities come from (01,02) = (4,0) for instance, and 
t~^ dampings are predicted. These discussions are summarized in Table [5l 

We explain the whole cancellation in {A4) (t) for (oi, 02) = (0, 0). The exact solution 
of (^44) (t) is proportional to 

(A,) {t) = C j giJu J2) sin((l]i(Ji, J2) - ^2iJu J2))t)dJ. (74) 

Exchanging the dummy variables Ji and J2, and using invariance of g under this 
exchange, we have (^44) (t) = 0. From the view point of singularity analysis, using 



notations of [Appendix D| and [Appendix E[ the invariance of g with respect to the 



exchange between Ji and J2 implies that the numerator of the integrand in ( ]E.3|) is 
of the type M"ln|M| with even a, since u = — J|. The evenness of a gives the 
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(01,02) 


(0,0) 


(1,0) 


(2,0) 


(^1) it) 




r2(N) 




{A2) it) 




r3(N) 


r3(c) 


{A3) it) 




t-3(N) 


r3(c) 


(A,) it) 


c 


r2(N) 


r2 



Table 5. Dampings and cancellations for the critical singularity. "C" means that 
there is a cancellation at leading order, and "N" that the leading order singularity 
does not exist. Dampings associated with the next leading singularities are observed 
when "C" or "N" appears. In {A4) (t) for (^1,02) = (OjO); the singularities at any 
orders are cancelled due to a special symmetry of g{m^ J). 

10 
1 

0.1 

0,01 

^ (1.001 
— 0.0001 
le-05 
le-06 
le-07 




0.1 
0.01 
0.001 
0.0001 
le-05 
le-06 
le-07 



ai = 1 










ai=2 J 


(b) 


i 



Figure 4. Dampings for the critical singularity. Obscrvables arc (a) and (b) 
A4. The parameter 02 is fixed as 0. The curves for ai = 1 and 2 in (b) are almost 
superposed. The curve of oi = in (b) is out of range due to a special cancellation. 
Except for this special cancellation, Ai and A2 give qualitatively similar panels to A4 
and A3 respectively. 



same sign for the two modes m and —m in the saddle singularity, hence they cancel 
each other in the sine observable ^44: this is of course compatible with the exact result 
{A,){t)=0. 

6.4- Toy model 

We set the frequencies as 

n^(j) = -j^-^- 2J2, n^iJ) = -2Ji + 2J2. (75) 

This toy model (!75l) has three types of singularities: 

• Vertex singularity at ( J^, J|) = (0, 0) for the modes (1, —1) and (—1, 1). 

• Tangent singularity at (J^, J2) = (1,0) for the modes (1, —1) and (—1, 1). 

• Line singularity on the J2 axis for the modes (1,1) and (—1,-1). 

No critical singularity and no singularity at infinity appear. We use distribution ( l66l) . 
with hi = h2 = ai = a2 = 0, that is 

ig{m, J) = e-(-^i+-^2)/2_ ^^g^ 
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V CI 


'T^Q 1^ O" 11 + 


T.iiic*! It^dfimo'i 


T iin 1 Gif^rrin H 1 




I 




















(0,0) 


(1,0) 


J2 axis 


J2 axis 


Ai = cos(^i + 62) 


N/A 


N/A 


c 


D 


A2 = sin(^i + 62) 


N/A 


N/A 


D 


C 


A, = cos{e^ - 62) 


E 


D 


N/A 


N/A 


= sin(6'i - 62) 


C 


D 


N/A 


N/A 



Table 6. Theoretical prediction of the asymptotic damping for the Ai observables 
in the toy model ((75|). "D" indicates the dominant singularity, "E" an existing 
but subdominant singularity, and "C" a cancelled singularity. "N/A" tells that the 
corresponding singularity cannot be captured by the observable. An oscillation is 
predicted for the tangent singularity only. 



6.4.1. Cancellation and higher-order contribution Three different types of singularities 
appear for the A^ observables, the strongest of which is a line singularity, see Table |6l 
Concentrating on Ai and A2, we show an example of cancellation for this line singularity. 

The relative sign between ^9(0;; m) and ^9(0;; — m) is (— l)-*^"^"^ for the line singularity. 
Thus, the leading singularity corresponding to ai = is cancelled for Ai but not for A2. 

Higher-order singularities may come from the expansion of g around Ji = and 
correspond to ai > 1. However, since g is even in Ji, the same cancellation will take 
place at all orders. Consequently, there is no higher-order contribution from g. 

Nevertheless, a higher-order contribution comes from m ■ f2(J). For the mode 
m = (1, 1), the function ip{uj; m) becomes 



(^(a;;(l,l)) 



1 

4i 



d.Ji 



c 



-3 Ji - Jf/2 - 00 

keeping the leading order of g{m, J), which is constant 
u = 3Ji + Ji/2, we have 

-C 



dJ' 



4i 



dJi 



-3 Ji - .Jl/2 - oj 
Changing variable from Ji to 



(77) 



I2i 
-C 



(1 + 2n/9)-V2 



I2i L 



u + UJ 



du 



-C 



12i 



1 - u/9 + O(n^) 
u + u 



In \uj\ In \uj\ 

I I 9 I I 



du 



A similar computation gives 
(/.(a;; (-1,-1)) 



-C 



12i 



'ill 1 I I 
m \uj\ m + 



(79) 



Coming back to the expression fjTOj) . the leading order ln|a;| is cancelled for the 
observable Ai = cos(6'i + 62). However, the second leading order a;ln|a;| survives, 
and gives a asymptotic damping. 

Similarly, we can find another cancellation for the observable A4 = sin(6'i — 62) on 
the vertex singularity. From the above discussions, we build Table [6] which summarizes 
existing, dominant and cancelled singularities for the four observables. 
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6.4-2. Numerical computations of the exact solution On figure \5[ we compare the 
theoretical prediction of Table [6] with the numerically computed exact temporal 
evolution of the four observables ( 168|) . We examine both amplitude's damping and 
frequency of (Aj) (t), whose exact solution is 



(A,) (t) 



TT 



e-iJ?+Jl)/^A,{ni{J)t,Q2{J)t)dJ,dJ2 



(80) 



by ( l65l) . Temporal evolutions of the four observables are reported in figure [5l and the 
agreements in damping rates are very good for all observables. Observables Ai and 
A2 have no tangent singularity, and hence the prediction for the frequencies of (Ai) (t) 
and (^42) (t) is zero. These zero frequencies are observed in figures [5]( a) and (b). On 
the other hand, the observables ^3 and ^4 are dominated by the tangent singularity 
at (J*, J2) = (1, 0), and the prediction for the frequencies is 0) — ^2(1, 0)| = 0.5. 

Looking at the temporal evolutions in linear scale in figure |6l we see that the predicted 
frequency is in good agreement with the numerics. 




T 



n.oi 
n.ooi 
o.noni 

le-05 



(b) 




1 10 ino 

t 



Figure 5. Temporal evolutions of the four observables: (a) {Ai){t) with Ai = 

cos(0i+02). (b)-(A2)(i)withA2 = sin(0i+02). (c) |(A3)(i)| with A3 = cos(0i-02). 
(d) |(A4)(t)| with A4 = sin(6'i— 02)- All panels have the same scales. Straight lines show 
the theoretically predicted algebraic dampings, which are respectively t~'^ ,t~^ ,t~^'^ 
and t~^-^. {Aj){t) is computed by using (|80| with cutoff at Ji = J2 = 10- 
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t 
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80 100 




Figure 6. Temporal evolutions of (a) {A3) and (b) {A4) in linear scale. Red plus 
points are numerically obtained, and green lines represent 4i~^-^ cos(0.5(t — 4.8)) and 
4t~^-^ cos(0.5(t — 8)) in (a) and (b) respectively. 



7. Spherically symmetric stellar system 

7.1. Isochrone model 

We present the isochrone model, which is a simple example of a spherically symmetric 
stellar system. The potential of this model is 

GM , . 



6 + + r2' 

where G is the gravitational constant, M is the total mass of the system, and 6 is a 
parameter. A stationary state satisfying the self-consistent condition for the potential 
is known to be jllj : 



/o 



v^(27r)3(GMfe)3/2 [2(1 - E)\^ 



+64E^ + 3(16E^ + 28E - 9) 



27 - 66^ + 320^2 - 240^3 



sm 



12) 



where E = —Eb/GM and E is energy. 

Let {r,^,ip) be the 3-dimensional polar coordinate, and {Jr,J^,J^) the conjugate 
angular momenta. We introduce new actions (Ji, J2, J3): 

Ji = Jip^ J2 = J-d + |^!/,|, ^3 = Jt- (83) 

We denote the conjugate angle variables by (^1,^2)^3)- The Hamiltonian of this model 
is then 

{GMf 



HiJ) 



84) 



The Hamiltonian is defined on the domain {Ji, J2, J3) G [— J2, </2] x [0, 00) x [0, 00). 
We remark that the Hamiltonian depends on J2 and J3, but not on Ji. We hence 
omit the modes corresponding to Ji, and write m = (777,2, ^3) for the mode vector 
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and f2 = (^2)^3) for the frequency vector. Thus, the isochrone model fits in the 2- 
dimensional analysis developed in this paper. Note that for consistency with standard 
notations, the actions J2 and J3 in the isochrone model correspond to Ji and J2 in the 
general theory of section [31 Frequencies are given by 

^3(J) = —{J) = ^ ' (85) 



.h + \{.h + ^JlTAGMh) 



njj) = ^(j) = -li+ ^ '^^ 1 risiJ). (86) 

7.2. Advection equation corresponding to isochrone model 

Using the frequencies Q2 ( l86l) and ^3 ( l85l) . we consider the advection equation 

This advection equation is obtained by omitting the potential perturbation in the 
linearized Vlasov equation. We take as initial state: 

f{e, j,t = o) = fo{J) + fi{e, j,t = o), (88) 

where fo{J) is the stationary state fl82|) and the perturbation fi is set as 

/i(6>, J,t = 0) = afo{J) cos(n2^2) cos(n3^3), |a| < 1. (89) 

Note that this perturbation is also independent of 61 and Ji. The Fourier transform 
ig{m, J) of fi{6, J,t = 0) is then 

zg{m, J) = I = (90) 

1 otherwise. 

Let us consider the temporal evolution of the expected value of an observable 
^(^2,^3)- As discussed in section 15.21 the Laplace transform of the expected value 
is given by ( l63i) . but the integrals should be performed over ^1,^2,^3,-^17-^2 and J3. 
Remembering that A, fi and f2 do not depend on 61 and Ji, the expected value may 
be written as a linear combination of functions of the form 

POO POO ^ j-J2 

ip{u]m) = / dJs / dJ2 / h{m, J)dJi. (91) 



^0 



m ■ Q,{J) — cu J^j^ 

Note that if we were studying the true hnearized Vlasov equation, the function h in the 
above equation would be replaced by m-'V jfo{J)ci{m, J)ck{m, J) or g{m, J)ci{m, J), 
according to the definitions of Fik{uj) and Gi{u) in equations (|25|1 and fl26|) respectively. 



Since hdJi vanishes for J2 = 0, a J2 factor can be factorized, and we can write 
J^^j^h{m, J)dJi = J2h{m, J). Thus, the functions fl9T|) fit in the framework of the 
abstract setting (1501) . and we apply the theory developed in sections [3] and HI with 

u{J) = J2him,J), fi{J) =m-n{J). (92) 

The J2 factor in i^( J) will be of importance to determine the singularity associated with 
vertex, line and infinite singularities. 
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Singularity 


Damping 


Frequency 


Condition for mode 


Vertex 




f 1712 1 ^ \ V GM 

I 2 +^5) 24fe3/2 


None 


Tangent 


^-1.5 


m2n2iJ*2,0)+mMJ*,0) 


—1 < 7723/7712 < —1/3 


Line 


r2 





77i2 = —27773 


Infinity 







None 



Table 7. Dampings and frequencies associated with singularities in the isochrone 
modeL This table shows the leading singularity for each type. (7712,7713) represents a 
mode. J| is the solution to ([94|) . Modes satisfying m2 + 27713 = or TO2 + 8777,3 = 
do not yield a generic vertex singularity, but a line singularity and a special vertex 
singularity respectively. Sec text for the leading line singularity. 



7.3. Singularities in isochrone model and theoretical prediction 
We can now list the singularities in the isochrone model: 

• Vertex singularity at (^2) -^s) — (OiO); except for the modes m satisfying m2 + 
2r?73 = or r?72 + 3r?73 = 0. The former case results in a line singularity, and the 
latter case gives a special vertex singularity with d{m ■ 17) / 9 J2 = at the origin. 

• Tangent singularity at (J25O) ^^e modes m satisfying —1 < 7773/7712 < —1/3. 
For such a mode, the singular point is given by the solution to the equation 

A m-n = 0, (93) 

which reads 



1 + \ 2^M^^7n3 

y/Ji + AGMbJ 3Jl + 4GMb m^' ^ ^ 

We remark that the left-hand-side of ( |94l) is a decreasing function of J2, and its range 
is (-1,-1/3) for J2 G (0,00). Thus 7773/7712 must be in the interval (-1,-1/3). 
The origin, namely J2 = 0, results in the special vertex singularity for the modes 
satisfying r?72 + 37773 = 0. 

• Line singularity on J3 axis for the modes m satisfying 7772 = —27773. 

• A singularity at infinity appears for any mode m. Estimations of z^(J^) and /i(J) 
are 

HJ) ~ (^2 + ^3)^', f^{J) ~ {J2 + J3)"', (95) 

since f22 ~ ^^3 ~ (J2 + -'3)"'^, fo ~ E^^'^, -E ~ (J2 + -^3)"^, and taking into account 
the J2 factor for u appearing in fl92l) . 

No critical singularity appears in the isochrone model. These singularities are arranged 
in Table [7] with asymptotic dampings and frequencies. 

The leading dampings are t~^, t~^'^, and for the vertex, tangent, line and 
infinity singularities respectively. We remark that the leading order of the function u is 
~ J2 around J2 = due to the J2 factor of u in (!92l) . thus the leading damping rates 
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for the vertex and line singularities respectively are not and t~^, but and . 
The dominant, i.e. slowest, decay is t~^/^ with zero frequency. 

In the next subsection we will observe the temporal evolution of the expected value 
of the observable A{6) = sin(n26'2 + n392) with respect to the exact solution 

h{e, J, t) = fi{e- n{j)t, j,t = o). (96) 

The mode {n2,n^) corresponds to the mode of perturbation ( l89ll . The expected value 
is: 

{A), (t) = j A{e)h{e, J, t)dedj. (97) 

For our interests, the prefactor of {A)^ (t) is not crucial; therefore we redefine (A)-^ (t) 



as 



{A)At)= / / J2fo{J)Min2^2iJ)+ns^s{J))t)dJ2dJ^ (98) 
Jo Jo 

in the following. For this sine observable, the leading vertex, tangent and infinity 
singularities survive, but the leading line singularity is expected to be cancelled looking 
at Table [21 since the corresponding exponents are ai = 1 and a2 = 0. From the line 
singularity, therefore, the second leading damping should arise instead of the leading 
damping t~^. However, since the line singularity is never dominant, this cancellation 
does not affect the following discussions. 



7.4- Numerical check 

We choose three modes: {n2,n3) = (1,1), (2,-1) and (3,-2). The mode (1,1) has 
vertex and infinity singularities, (2, —1) has all four singularities, and (3, —2) has all 
four except the line singularity. To perform the numerical integration of (p8l) . we set 
the parameters asG' = M = 6 = l. The exact solutions ( |98l) are shown in figure [7l 

The mode (1,1) shows t~'^^^ damping without oscillation as the theory predicted. 
This mode also has a vertex singularity, which gives a non-zero frequency, but no 
oscillation is observed in the asymptotic time region: the vertex singularity contribution 
is expected to damp as t~^, which may be too fast to be visible at large times. 

The modes (2, —1) and (3, —2) asymptotically damp as as predicted, but they 
have oscillations. At variance with the mode (1, 1), both these modes have a tangent 
singularity, which gives a rather slow damping with non-zero frequency. To confirm 
that there is a contribution from the tangent singularity, we show the power spectra of 
(A)^ (t) in figure [HI The peaks of the power spectra are in good agreement with the 
theoretically predicted frequency |m2n2(J2,0) + m3r23(J2 , 0)|, where J2 is the solution 

to dnn). 

As mentioned above, to study the true Vlasov equation, we should also estimate 
the Cfc(m, J) functions, which depend on the biorthogonal functions Uk{q)- 
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Figure 7. Damping in the isochrone model for the observable A = sin(n26'2 + n^Os). 
("-2, "-3) = (111) (red plus), (2,-1) (green cross) and (3,-2) (blue star). The purple 
guide line is proportional to t^^^^. {A){t) is computed by using (|98| with cutoff at 
J2 = J3 = 20. 
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Figure 8. Power spectra of {A)i{t) for the modes (2,-1) (green cross points) and 
(3,-2) (blue star points). The two vertical lines indicate to = 0.1185 and 0.0509, 
which are predicted theoretically for the modes (2,-1) and (3,-2) respectively as 
|m2fi2(>/2 ' 0) + "T.3r23( J2 , 0)1 with J2 the solution to (|M|) . The power spectra are 
computed from time series with time slice 1 in the time interval [77, 1100] to avoid an 
early transient time region. 



8. Conclusion 

When a stable inhomogeneous stationary state of a Vlasov equation is perturbed, 
the asymptotic damping of the perturbation is algebraic, no matter how regular 
the stationary state and perturbation are. The damping rate and frequency are 
controlled by the singularities on the real axis of the Fourier-Laplace transform of the 
perturbation. In this paper, in systems with two spatial dimensions, we have classified 
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^ K,>0 



K|=0 



Figure Al. The boundary of D is regular close to point A, its intersection with the 
curve /x(Ji, J2) — xq, that is Ki — 



these singularities into vertex, tangent, critical, line and infinity singularities, and have 
given damping rate and frequency in each case. This classification is also valid for 
three dimensional spherically symmetric stationary states, since they reduce to the two- 
dimensional case thanks to symmetry. The resulting picture is much richer than in 
one spatial dimension [22]. We have illustrated the theory on a toy model, and on 
an advection equation associated with the isochrone model, a simple model for self- 
gravitating systems: tests in these models show that the singularity analysis captures 
very well the observed damping. 

The main goal of this theory is a better understanding of the asymptotic relaxation 
of self gravitating systems. To achieve this goal, a detailed study of the classical 
expansion of a perturbation on a biorthogonal basis should be coupled to the singularity 
analysis we have presented here. 
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Appendix A. Analytic continuation 

The starting point is the following formula for (p{z), valid for lm{z) > 0: 



We want to understand the regularity of ip{z) close to 2; = G M. 

Assuming that Xq is not a critical value for the function /x (condition i), the equation 
f^iJi, J2) = Xq defines one or several curves on the (Ji, J2) plane. We treat here the case 
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with one curve, the argument carries over easily to the case of several curves. 

It is then possible to use a new set of variables Ki, K2 fulfilling the two conditions: 

a) Ki = /i(Ji, J2) - xo 

b) K2 is such that the Jacobian dKidK2/dJidJ2 never vanishes. We do not need to 
further specify K2. 

Then, with a suitable definition of u: 



If the boundary of the domain D is regular around its intersection with the curve 
Ki = (condition ii) (see figure \KT\ point A), then K2^^\Ki) and K2"^\Ki) are regular 
functions of Ki , and we are left with the single integral 



where /i is a regular function. The usual Landau continuation argument then ensures 
that (p{z) is regular at 2; = Xq, unless k["^^ = or k[^^ = (condition iii). 

Thus, under conditions i), ii) and iii), ^p{z) is not singular a.t z = xo- Breaking 
one of these conditions yields a critical (section [375]) . vertex (section 13731) and tangent 
(section 13.41) singularity respectively. 

Appendix B. Vertex singularity 

Let us compute the leading singularity of 



The origin {u,v) = (0,0) is on the boundary of the domain U, and is a singular point. 
We assume that the boundary of the domain U is locally approximated by two lines 
around the origin, and that the line u = does not coincide with one of such boundary 
lines. We can classify the shapes of the domain in three cases: (i) the two lines are in 
the u > half plane; (ii) the two lines are in the m < half plane; (iii) one line is in the 
u > half plane, and the other in the u < half plane. Case (ii) is made identical to 
case (i) by changing the sign of u, z and C if necessary. Thus we consider case (i) and 
case (iii). Remembering that the line u = does not coincide with the boundary of f/, 
we can further classify V into the following four types: 




Integrating over K2 first: 






(B.l) 



{(M,f) I < M < e, f3u < V < au}. 



(B.2) 



U2 



{{u,v) I < M < e, au<v<6oT6<v< /3u} 
Li{{u,v)\ - e < u < 0, -5<v<5}. 



(B.3) 
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Figure Bl. Four types of domain U for the vertex singularity. 




u 



U3 = {{u,v) I < n < e, au < V < 6} 

U {(m, I - e < u < 0, I3u < V < 6}. 



(B.4) 



{{u,v) I < M < e, 6 < V < au} 
U{{u,v) I - e < n < 0, 6 <v < f3u}. (B.5) 

The domains Ui and U2 correspond to case (i), and U3 and U4 to case (iii). 6, e, a and 
/3 are constants. See figure [BTI for schematic pictures of the domains. We first note that 
the integral over the domain 

f/rcgular = {{u,v) \ - € < U < €, aU < V < 6} (B.6) 

does not give any singularity for the function (IB. II) . This fact is immediately obtained 
by performing the integral over v. 



u^du 



V dv 



u 



- [au] 



i+i 



du 



(B.7) 



U-{z- Xq) U-{z-Xq) l + l 

and using Appendix E[ Adding or subtracting such a regular domain, domains U2,U3 
and t/4 can be reduced to f/i. Thus, what we have to consider are integrals of the form: 



u^du 



u — {z — Xo) J 



i/dv 



a 



l + l 



u 



k+l+l 



U-{z-Xq) 



du [B. 



Using [Appendix E| again, the leading singularity of = limy^o+ '^{x + iy) at xq is 

(j)i:^^{x) = ci(x - xo)''+'+Mn |x - xol + C2(x - x^f+^^^H^x - Xq) (B.9) 
except in the non-generic situation (3 = —a, I odd. 



Appendix C. Tangent singularity 

Changing variables to u = Jf + J2,v = Ji, the leading singular part of the function (p{z) 
is: 

du 



^{z) = C 



u-{z-xo}J^^ 
"° u'^^-'^du 



v^'iu-v'Y^dv 



u-[z-Xo)j_^ 



(C.l) 
(C.2) 



On algebraic damping close to inhomogeneous Vlasov equilibria in multi- dimensional spacesZl 
The integral over v is computed as 



V 



ai+2k 



dv = = (1 - (_l)l+«i+2fe)^(l+ai+2fc)/2_ 



^ 1 + ai + 2k 



Thus we have 



I _ r_iy+ai+2k ruo ^l/2+ai/2+a2 

The prefactor 1 — (_i)i+ai+2fc = i _ (^_\Y+ai vanishes for any k when ai is odd, 
thus we may assume that Oi is even. We are left with integrals as in flE.2p . with 
1/2 + ai/2 + 02 ^ Thus the leading singularity is 

=C'i(x-xo)'/'+"^/'+''^if(x-a;o)+C72(xo-x)^/2+"^/'+'^2i/(xo-x).(C.3) 
Appendix D. Critical singularity 

Calling Ai and A2 the two eigenvalues of the Hessian of yU at the critical point, we expand 
/i( J) after an appropriate change of basis 

/i(J) ~ xo + Ai(Ji - Jlf + A2(J2 - j;?. (D.l) 

We have to consider two cases: i) the critical point is a local extremum; ii) the critical 
point is a saddle. 

Case i), A1A2 > : Without loss of generality, we assume Ai > , A2 > 0. Using the 
leading order approximation ( IMl) for u and performing the usual shifting and scaling of 
variables, we obtain the leading singular part of ^{z) as 

r ja2 

2tt ^l+ai+a2 







C / cos"^ a sin"2 arfa / — -dr (D.3) 

Jo H - (2 - xq) 



where (r, a) are the polar coordinates on the ( Ji, J2) plane. A change of variable u = 
reduces the integral over r to the form (IE.2p and we conclude that the leading singularity 



IS 

' Ci{x - XoY''''+''^'^/^\n\x - Xo\ 

+C2ix - Xo)("^+"2)/2i/(x - Xo) fli, aa : even (D.4) 
otherwise. 

Case ii), X1X2 < : 

Without loss of generality, we may assume, shifting and scaling Ji and J2, that 
Ai = 1 , A2 = — 1, and the saddle is at {Ji, J2) = (0, 0). We need then to study 

/Cl /•C2 jai ja2 
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The function (f{z) is zero except for even Oi and 02- We perform the change of variables 
ji = Ji + J2 , j2 = Ji- J2 ; dji dj2 = 2 dJi dJ2. (D.6) 

Forgetting a constant factor and modifying the integration domain (the important thing 
is to integrate over a neighborhood of (0, 0)), we are led to study 

<f{z) = / / — 7 ^ — djidj2 

"1 °2 ^ci 1-02 -k+l ■ai+a2-k-l 

Jl J2 



k=0 1=0 

After a further change of variables 



^ ^ J-ci J-C2 J1J2 -{Z- Xo) 



u = jij2 , V = j2 ; dji dj2 = -^du dv (D.8) 

\v\ 



we obtain 



ai a2 
k=0 1=0 

(/•C1C2 I-C2 pO j-C2 \ ^k+l^ai+a2-2{k+l) 

I du dv+ du I dv] -— -. (D.9) 

Jo Ju/c^ J-c,C2 J-u/c^ J {u-{z-Xo))\v\ 

Remembering that ai + 02 is even, we perform the integration over v: 

"1 "2 ^ClC2 ^k+l 



J. n 1 n J —( 



k=0 1=0 >'-ciC2 

i(c] - (n/ci)^), 7 = ai + a2 - 2(A; + /) 7^ 
In IC2I — In |n/ci|, 7 = ai + 02 — 2{k + /) = 0. 



(D.IO) 



The terms with 2{k + I) 7^ ai + 02 are not singular, so we are left with 



I u — z 



k+l={ax+a2)/2 



According to Appendix E again, the singularity of is then 

=Ci(x-Xo)("^+"^^/Mn|x-Xo|+C2(x-Xo)("^+''^^/^i/(x-Xo)(D.12) 

for ai and 02 even. 

Appendix E. Singularity computation 

We first show that the function 

(hix) = lim / -du, a G Z . (E-1) 

y^o+J^,u-{x + iy) 

does not have any singularity. Then we compute singularities of the functions 

= lim / -du, a G M\{---,-3, -2, -1}, (E.2) 

y^0+ Jq u- {x + ly) 
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and 

mix) = lim / -du, a E Z 4. . (E.3) 

y^o+ J^^u-{x + iy) 

The constant c is assumed to be positive in the above three functions. All the 
singularities found in this paper reduce to a computation of the singularity close to 
a; = of one of the above two functions, (IE. 21) and (IE. 30 . We show in this appendix 
that singular parts of are 



ismg 



[x) = — In \x\ + Z7rx"if(x), a G Z_| 



(E.4) 



and 



= Cix"H{x) + C2{-xYH{-x) + i7rx"H{x), a G M \ Z (E.5) 
for the function ( 1E.2|) . and 

0«i^s(x) = Cx"H{x) + ZTTX^ln a G Z+ (E.6) 
for the function (lE.Sp . 



Appendix E.l. No singularity of the function liE. 1\) 
We start from the function (IE. II) 

{u — x)u' 



(t){x) 



lim 



-du + 



:du 



(E.7) 



We denote the real and the imaginary parts of by (j)R{x) and respectively. 
The imaginary part can be computed using the change of variable s = {u — x) / y: 

(j)i{x) = 7rx''{H{x + c) - H{x-c)), (E.8) 

where H is the Heaviside step function. Thus 4'i{x) has no singularity around a; = 0. 
The real part is simply: 



fx) = PV 



PV 



u 



U — X 

X' 



-du 



/C G. 1 /i 
du + Ci / 
„ n — X ^-^ ' J_ 



\u — X] 



\a-l-\ 



du 



(E.9) 



where PV denotes the principal value. The sum on the right-hand-side is clearly regular, 
and the first term gives: 



PV 



X 



-du = x"(ln |c — x| — In I — c — x|). 



u — X 

so that no singularity appears around x = 0. 



(E.IO) 



On algebraic damping close to inhomogeneous Vlasov equilibria in multi- dimensional spaces'i'i 
Appendix E.2. Singularity of the function liE. S\) 



We study the function (lE.2p . which differs from the function ( lE.ll) in the lower bound 



of the integral. The imaginary part, (pi, is then: 

(f)i{x) = 7rx''{H{x) - H{x-c)), (E.ll) 
and hence the singularity of (pi around x = is 

0f = 7rx"i7(x). (E.12) 
The real part is also directly obtained using equation (lE.lOp : 



= x"(ln|c-x| -ln| (E.13) 

for a G Z4.. Thus, the singularity of (pR around x = is 

(P'P{x) = -x" In |a;|, a G Z+. (E.14) 

Combining real and imaginary part gives (IE.4I) . 

The final step for the function (1E.2I) is to investigate the real part for a G M \ Z. 
We divide (Pr{x) in two parts: 

(Pr(x) = lim / du+ / du . (E.15) 

e-^O U-X 7,^1+, U-X J 

This division is just the definition of the principal value for x > 0, and is also valid for 
X < 0. We then use the expansions, for m < |x| and m > |x| respectively: 
1 -.00 , -, ^00 , 

1 -lE(-) ; — = -E(-) ■ (Ew) 

■ ^-^ \x/ u — x u^-^\u/ 



U — X X „ „„„ 

A;=0 fc=0 



Substituting the above expressions into flE.lSp and remembering a G M \ Z, we have 

1^ I l^k=0\a+k+l^ a-k ) ^ 2^k=0-^ a-k ^ U. 

Note that all series converge. The second series for x > and x < exactly coincide, 
hence they do not contribute any singularity. The singularity comes from the first series 
and is: 



x) = Cix°if(x) + C2(-x)"//(-x), aGM\Z. (E.18) 



Consequently, we have proved flE.5l) for the function flE.2p 



Appendix E. 3. Singularity of the function liE. gj) 

The imaginary part of the function ( 1E.3|) . denoted by (pi{x), is 

(piix) = 7rx"ln|x|(i7(x + c) -i7(x-c)), (E.19) 
and the singularity around x = is 

= 7rx°ln|x|. (E.20) 
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To compute the real part, we rewrite the integrand as done in (1E.9P : 



'x) = PV 



X 



u — X 



^o.Ci{u-xY-^-^x^ 



1=0 



In 



We remark that the integral 

n" In \u\du {n G Z+) 



(E.21) 



(E.22) 



can be performed in the sense of improper integral, and converges to a finite value. 
Thus, concentrating on the singularity, the real part is reduced to 



In \u\ 

(x) =x'' PV I ^^du 
_^u-x 

hm ( / + 



\x\+e 



In lul 



-du. 



\x\+e^ 



U — X 



(E.23) 



Using expansion (lE.16p . we obtain 



1^1-^ ln|M| , ^ -2 

du = y 

.|^|+, u-x ^ 



and 



X 

k=0 

-1^1-^ r \ ln|M| 

J\x\+el U-X 



(kl-e) 



2k+l 



2k + 1 



In I Ixl — el 



\x\ — e 



v2fc+l 



(2A; + 1) 



(E.24) 



du 



k=0 



X 



2k+i ( +e) 



-(2fc+l) 



2k + I 



In I Ixl 



{\x\+e) 



-(2fe+l) 



{2k + I] 



(E.25) 



where we have omitted regular functions in (IE. 251) . The logarithmic terms are cancelled 
by adding (IE.24p and (1E.25P and taking the limit e — > 0; hence the singularity of (pnix) 
is: 



= Cx'^Hix). 



(E.26) 



The imaginary part ( 1E.20p and the real part (]E.26p prove (1E.6P for the function ( lE.Sp . 



Appendix E.4- Relative sign between modes m and —m 

Coming back to the functions F's and G's, we are interested in singularities of the 
function 



ip{z;m) 



-dJ. 



(E.27) 



m ■ fl{J) — z 

To discuss a possible cancellation between the modes m and — m, we show a relation 
between the above function and 



ip{z\ -m) 



-dJ. 



(E.28) 



-m ■ r2( J) — z 

We note that the numerator J) may also depend on m; in this case one would have to 
discuss the sign of the numerator separately; we ignore this dependence in this section. 
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Remembering that the singularity of the function (1E.27P results in one of the two 
types of functions flE.2p and flE.Sp . we consider the relation between 

and 

<^"xn(^) = - lim f _clu (E.30) 

where h{u) is (a G Z+ U (M \ Z) or In \u\ (a G Z+), and = m ■ ri( J*) e M with 
a special point J = J* . (pZ^^ is obtained by changing the signs of the prefactor, x and 
y in 0+ , thus the singular points of 4>toi^) ^-xoi^) ^^'^ respectively. If 

Xo 7^ 0, no cancellation occurs in general due to the different singular points. If Xq = 0, a 
cancellation may occur depending on the relative sign of the singularities for the modes 
m and — m. The change of sign of y implies the change of sign of the imaginary part, 
and we use H{—x) = 1 — H{x). Hence the relation between the singularities of 
and 00 (x) is: 

00 "'"'(^) = {-lf^''<Pr'^\x\ a G Z+ (E.31) 
for the singularity (IE. 40 . and 

00 •^'"'(^) = (-l)>o'"'''"'(^), « e Z+ (E.32) 
for the singularity (IE. 60 . The relations (lE.Sip and (1E.32P imply that, for instance, 
singularities of ip{z\ m) + ip{z; —m) and of ip{z] m) — ip{z; —m) cancel respectively if a 
is or positive even and xq = 0. There is no simple relation such as (lE.Sip for singularity 
(lE.Sp . thus no cancellation is expected in general. 
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